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Summary 

A transonic unsteady aerodynamic and aeroelasticity code called CAP-TSD has been developed for 
application to realistic aircraft configurations. The name CAP-TSD is an acronym for Computational 
Aeroelasticity Erogram - Xransonic Small Cisturbance. The code permits the calculation of steady 
and unsteady flows about complete aircraft configurations for aeroelastic analysis in the flutter 
critical transonic speed range. The CAP-TSD code uses a time-accurate approximate factorization 
(AF) algorithm for solution of the unsteady transonic small-disturbance potential equation. The 
paper gives an overview of the CAP-TSD code development effort and reports on recent algorithm 
modifications. The algorithm modifications include: an Engquist-Osher (E-O) type-dependent 

switch to treat regions of supersonic flow, extension of the E-0 switch for second-order spatial 
accuracy, nonisentropic effects to treat strong-shock cases, nonreflecting far field boundary 
conditions for unsteady applications, and several modifications to accelerate convergence to steady- 
state. Calculations are presented for several configurations including the General Dynamics one- 
ninth scale F-16C aircraft model to evaluate the algorithm modifications. The modifications have 
significantly improved the stability of the AF algorithm and hence the reliability of the CAP-TSD 
code in general. Calculations are also presented from a flutter analysis of a 45° sweptback wing 
which agree well with the experimental data. The paper presents descriptions of the CAP-TSD code 
and algorithm details along with results and comparisons which demonstrate the stability, accuracy, 
efficiency, and utility of CAP-TSD. 


Introduction 

Presently, considerable research is being conducted to develop finite-difference computer codes 
for calculating transonic unsteady aerodynamics for aeroelastic applications. 1 These computer codes 
are being developed to provide accurate methods of calculating unsteady airloads for the prediction of 
aeroelastic phenomena such as flutter and divergence. For example, the CAP-TSD2 unsteady 
transonic small-disturbance (TSD) code was recently developed for transonic aeroelastic analyses of 
complete aircraft configurations. The name CAP-TSD is an acronym for Computational 
Aeroelasticity Erogram - Xransonic Email disturbance. The new code permits the calculation of 
unsteady flows about complete aircraft for aeroelastic analysis in the flutter critical transonic speed 
range. The code can treat configurations with arbitrary combinations of lifting surfaces and bodies 
including canard, wing, tail, control surfaces, tip launchers, pylons, fuselage, stores, and nacelles. 
Steady and unsteady pressure comparisons were presented in Refs. 2 and 3 for numerous cases which 
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demonstrated the geometrical applicability of CAP-TSD. These calculated results were generally in 
good agreement with available experimental pressure data which validated CAP-TSD for multiple 
component applications with mutual aerodynamic interference effects. Preliminary aeroelastic 
applications of CAP-TSD were presented in Ref. 4 for a simple well-defined wing case. The case was 
selected as a first step toward performing aeroelastic analyses for complete aircraft configurations. 
The calculated flutter boundaries compared well with the experimental data for subsonic, transonic, 
and supersonic freestream Mach numbers, which gives confidence in CAP-TSD for aeroelastic 
prediction. 

The CAP-TSD code uses a time-accurate approximate factorization (AF) algorithm recently 
developed by Batina 5 for solution of the unsteady TSD equation. The AF algorithm involves a Newton 
linearization procedure coupled with an internal iteration technique. In Ref. 5, the algorithm was 
shown to be efficient for application to steady or unsteady transonic flow problems. It can provide 
accurate solutions in only several hundred time steps, yielding a significant computational cost 
savings when compared to alternative methods. For reasons of practicality and affordability, an 
efficient algorithm and a fast computer code are requirements for realistic aircraft applications. 

Recently, several algorithm modifications have been made which have improved the stability of 
the AF algorithm and the accuracy of the results. 6 - 7 These algorithm modifications include: (1) an 
Engquist-Osher (E-O) type-dependent switch to more accurately and efficiently treat regions of 
supersonic flow, (2) extension of the E-0 switch for second-order-accurate spatial differencing in 
supersonic regions to improve the accuracy of the results, (3) nonisentropic effects to more 
accurately treat cases with strong shocks, (4) nonreflecting far field boundary conditions for more 
accurate unsteady applications, and (5) several modifications which accelerate convergence to 
steady-state. The work has been a major research activity over the past two years within the 
Unsteady Aerodynamics Branch at NASA Langley Research Center. The purpose of the paper is to give 
an overview of the CAP-TSD code development effort and report on the recent algorithm changes and 
code improvements. The paper documents these developments and presents results which 
demonstrate the capability. 


Symbols 


bo reference length, cr/2 

c airfoil chord 

unsteady lift-curve slope 

Cr wing reference chord 

Cp pressure coefficient 

k reduced frequency, wcr/2U 

M freestream Mach number 

NSUP number of supersonic points 

q i( t ) generalized displacement for mode i 
R residual 

t time, nondimensionalized by freestream speed and wing reference chord 

U freestream speed 

a instantaneous angle of attack 

ao, ai mean angle of attack and amplitude of pitch oscillation, respectively 
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y ratio of specific heats 

At nondimensional time step 

t\ fractional semispan 

p ratio of wing mass to mass of air in the truncated cone that encloses the wing 

p freestream flow density 

<{> disturbance velocity potential 

to angular frequency 

coa natural frequency of the first torsion mode 

Subscripts 

t tail 

w wing 


Transonic Small-Di sturbance Equation 

The flow is assumed to be governed by the general frequency modified TSD potential equation 
which may be written as 

M 2 (4>, + 2V, = [(1 - M 2 )^ + F* 2 + G<$ x + (<t> y + H^Vy + «>A < 1 ) 


Several choices are available for the coefficients F, G, and H depending upon the assumptions used in 
deriving the TSD equation. For transonic applications, the coefficients are herein defined as 

F = -i(y+l)M 2 , G=i(y-3)M 2 , H = -(y-1)M 2 (2) 

The linear potential equation is recovered by simply setting F, G, and H equal to zero. 


Approximate Factorization Algorithm 

A time-accurate approximate factorization algorithm was developed 5 to solve the unsteady TSD 
equation. In this section, the AF algorithm is briefly described. 


General Description 

The AF algorithm consists of a Newton linearization procedure coupled with an internal iteration 
technique. For unsteady flow calculations, the solution procedure involves two steps. First, a time 
linearization step (described below) is performed to determine an estimate of the potential field. 
Second, internal iterations are performed to provide time-accurate modeling of the flow field. 
Specifically, the TSD equation is written in general form as 
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(3) 


R ((J) n+1 ) = 0 


where <(> n+1 represents the unknown potential field at time level (n+1). The solution to Eq. (3) is 
then given by the Newton linearization of Eq. (3) about 

* 

R(<t>) + (— X *A<() = 0 (4) 

a<i> ** 

In Eq. (4), <t>* is the currently available value of <t> n+1 and A<t> = «j>n+l - $*. During convergence of the 
iteration procedure, A<j> will approach zero so that the solution will be given by <t> n+1 = <|>\ In 
general, only one or two iterations are required to achieve acceptable convergence. For steady flow 
calculations, iterations are not used since time accuracy is not necessary when marching to steady 
state. 


Mathematical Formulation 

The AF algorithm is formulated by first approximating the time-derivative terms by second- 
order-accurate finite-difference formulae. The TSD equation is rewritten by substituting <t> = <t>‘ + 
a<j> and neglecting squares of derivatives of A<{> (which is equivalent to applying Eq. (4) term by 
term). The resulting equation is then rearranged and the left-hand side is approximately factored 
into a triple product of operators yielding 


L L L A<J> = - a R(<t>*, <|> n , <)> n 2 ) (5) 

s i s 

where the operators L4, U|, L; and residual R were derived and presented in Ref. 5. In Eq. (5) a is a 

relaxation parameter which is normally set equal to 1.0. To accelerate convergence to steady-state, 
the residual R may be over-relaxed using a > 1 . Equation (5) is solved using three sweeps through 
the grid by sequentially applying the operators L4, L^, L; as 


£, - sweep: 

A<() = - G R 

(6a) 

T| - sweep: 

L A6 =A6 
*n 

(6b) 

£ - sweep: 

11 -e- 
< 
11 

-©- 

< 

(6c) 


Further details of the algorithm development and solution procedure may be found in Ref. 5. 
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Flow-Tanaencv Boundary Conditions 


The flow tangency boundary conditions are imposed along the mean plane of the respective lifting 
surfaces and the wakes are assumed to be planar extensions from the trailing edges to the 
downstream boundary of the finite-difference grid. The numerical implementation of these 
conditions 2 allows for coplanar as well as non-coplanar combinations of horizontal (canard, wing, 
horizontal tail, launchers) and vertical (pylons, vertical tail) surfaces. Bodies such as the fuselage, 
stores, and nacelles are treated using simplified boundary conditions on a prismatic surface rather 
than on the true surface. 2 The method is consistent with the small-disturbance approximation and 
treats bodies with sufficient accuracy to obtain the correct global effect on the flow field without the 
use of special grids or complicated coordinate transformations. This type of modeling is similar to 
that of Boppe and Stern, 8 where good agreement was shown in comparison with experimental data for 
configurations with a fuselage and flow-through nacelles. 


Time-Linearization Step 

An initial estimate of the potentials at time level (n+1) is required to start the iteration 
process. This estimate is provided by performing a time-linearization calculation. The equations 
governing the time-linearization step are derived in a similar fashion as the equations for iteration. 
The only difference is that the equations are formulated by linearizing about time level (n) rather 
than the iterate level (*). 


Algorithm Improvements 


Enqquist-Osher Type-Dependent Switch 

Algorithms based on the TSD equation typically use central differencing in regions of subsonic 
flow and upwind differencing in regions of supersonic flow. This, of course, allows for the correct 
numerical description of the physical domain of dependence. The original CAP-TSD code of Ref. 2 

used the Murman 8 type-dependent switch to change the spatial differencing. The Murman switch, 
however, admits nonphysical expansion shocks as a part of the solution and has been shown to be less 
stable than monotone methods.io.n For example, unsteady results for a NACA 64A006 airfoil were 
presented in Ref. 11 which demonstrated an order of magnitude increase in time step using a 
monotone algorithm. Therefore, an Engquist-Osher (E-O) monotone switch, similar to that of Ref. 
10, has been incorporated within the AF algorithm of the CAP-TSD code. 6 The E-0 switch is based on 
sonic reference conditions and does not admit expansion shocks as part of the solution. Use of the E-0 
switch also generally increases computational efficiency because of the larger time steps which may 
be taken. 


Second-Order-Accurate Spatial Differencing 

Most TSD algorithms are only first-order-accurate spatially in regions of supersonic flow. This 
is due to the first-order upwind differencing that is typically used to treat these regions. Use of 
second-order upwind differencing has been shown to improve the accuracy of the solution while 
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retaining the numerical stability of the first-order method. 12 Consequently, the E-0 type dependent 
switch of the AF algorithm has been extended for second-order spatial accuracy in supersonic 
regions of the flow. 8 Comparisons of results obtained using first-order and second-order 
differencing, to be presented, demonstrate the improved accuracy of the second-order method. 


Entropy and Vorticity Effects 

A serious limitation of the potential flow codes, in general, is the inability to accurately treat 
cases with strong shock waves. For these cases, the isentropic potential formulation typically 
overpredicts the shock strength and locates the shock too far aft in comparison with experiment. In 
fact, it is fairly well known that potential theory predicts non-unique steady-state solutions 13 * 14 
for narrow ranges of Mach number and angle of attack. Simple modifications to potential theory, 
however, have been shown to eliminate the nonuniqueness problem and consequently provide 
solutions which more accurately simulate those computed using the Euler equations. 15 ’ 18 These 
modifications include the effects of shock generated entropy and they require only minor changes to 
existing computer codes. 

Rotational effects may also become important when strong shock waves are present in the flow. 
For example, vorticity is generated by shock waves due to the variation of entropy along the shock. 
Potential theory, of course, does not account for these effects because of the irrotationality 
assumption necessary for the existence of a velocity potential. For these cases the Euler equations 
are generally required to accurately model the flow. Recently, however, simple modifications to 
potential theory have been developed to model rotational effects. 1 9.20 These modifications involve a 
velocity decomposition originally suggested by Clebsch .21 in this model, the velocity vector is 
decomposed into a potential part and a rotational part. For most applications of interest to the 
aeroelastician, the rotational part occurs only in the region downstream of shocks. Therefore, the 
potential part can be obtained throughout most of the flow field using an existing potential flow code. 
The rotational part can then either be incorporated as a source term in the governing equation or as a 
modification to the fluxes, by making simple changes to the solution procedure. These changes 
consequently include the effects of shock generated vorticity as well as entropy and require 
relatively straightforward modifications to existing potential codes. 

Entropy and vorticity modifications to TSD theory have been implemented within the AF 
algorithm of the CAP-TSD code as described in Ref. 7. These modifications include: (1) an 
alternative streamwise flux in the TSD equation which was derived by an asymptotic expansion of the 
Euler equations, 22 (2) a modified velocity vector which is the sum of potential and rotational parts 
which in turn modifies the streamwise flux, and (3) the calculation and convection of entropy 
throughout the flow field. The modified code includes the effects of entropy and vorticity while 
retaining the relative simplicity and cost efficiency of the TSD formulation. Results obtained 
including these effects, to be presented, clearly demonstrate the improved accuracy of the AF 
algorithm. 


Nonreflectina Far Field Boundary Conditions 

For unsteady applications, the far field boundary conditions can have a significant influence on 
the accuracy of the solution. Steady-state boundary conditions are inadequate for unsteady 
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calculations, since disturbances reaching the boundaries are reflected back into the computational 
domain. These reflected disturbances can propagate into the near field and thus produce inaccurate 
results. One solution to this problem is to locate the grid boundaries far away to minimize the effect 
of the boundary conditions. This is generally not an acceptable remedy because of the higher 
computational cost which results from an increased number of grid points required to discretize a 
larger computational domain. The more appropriate solution is the use of nonreflecting far field 
boundary conditions which absorb most of the waves that are incident on the boundaries and 
consequently allow the use of smaller computational grids. 23 Nonreflecting boundary conditions 
similar to those of Whitlow 23 have been incorporated within the CAP-TSD code. 6 These boundary 
conditions are consistent with the AF solution procedure and are described in more detail in Ref. 6. 
Results obtained with and without the nonreflecting boundary conditions are presented which 
demonstrate their effectiveness. 


Steadv-State Convergence Acceleration 

Finally, several algorithm changes have been made to accelerate convergence to steady-state. 6 
Besides the E-0 switch, these changes include: (1) deletion of the time-dependent terms from the 
residual of the AF algorithm, (2) deletion of all of the time-derivatives of the TSD equation, and (3) 
over-relaxation of the residual. The effects of these modifications on the steady-state convergence 
are demonstrated in the results presented herein. 


Aeroelastic Solution 

In this section the aeroelastic computational procedures are described including the equations of 
motion and the time-marching solution. 


Equations of Motion 

The aeroelastic equations of motion are based on a right-hand orthogonal coordinate system with 

the x-direction defined as positive downstream and the z-direction positive upward. 4 The 
presentation herein is limited to the case of an isolated wing with motion in the z-direction from an 
undisturbed position in the z = 0 plane. The equations of motion may be written as 

Mq + Cq + Kq = Q , 7) 


where q is a vector of generalized displacements, M is the generalized mass matrix, C is the damping 
matrix, and K is the stiffness matrix. In Eq. (7) Q is the vector of generalized aerodynamic forces 
defined by 



pU 2 2 
2 


J 


h. 

1 


S 


Ap dS 
pU 2 /2 c 2 


( 8 ) 
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where Ap is the lifting pressure which is weighted by the mode hj. The equations of motion (Eq. (7) 
are solved by time-integration by first rewriting Eq. (7) as 

q = - M 1 Kq - M 1 Cq + M 1 Q (9) 


Time-Marching Solution 

The aeroelastic solution procedure for integrating Eq. (9) is similar to that described by 
Edwards, et al . 24 Reference 24 describes for a two-dimensional, two-degree-of-freedom system an 
aeroelastic solution in terms of a state equation formulation. By a parallel formulation, a linear 
state equation is developed from Eq. (9) which is solved numerically using the modified state- 
transition matrix integrator of Ref. 24. This integrator was shown to be superior to six alternative 
structural integration algorithms. 2 ® 

For aeroelastic analysis, two steps are generally required in performing the calculations. In the 
first step, the steady-state flow field is calculated to account for wing thickness, camber, and mean 
angle of attack thus providing the starting flow field for the aeroelastic analysis. The second step is 
to prescribe an initial disturbance to begin the structural integration. Disturbance velocities in one 
or more modes, rather than displacements, have been found to be distinctly superior in avoiding 
nonphysical, strictly numerical transients and their possible associated instabilities. In 
determining a flutter point, the freestream Mach number M and the associated freestream speed U 
are usually held fixed. A judiciously chosen value of the dynamic pressure p(J 2 /2 is used to compute 
the free decay transients. These resulting transients of the generalized coordinates are analyzed for 
their content of damped or growing sine-waves, with the rates of growth or decay indicating whether 
the dynamic pressure is above or below the flutter value. This analysis then indicates whether to 
increase or decrease the value of dynamic pressure in subsequent runs to determine a neutrally 
stable result. 


CAP-TSD Code 

The AF algorithm is the basis of the CAP-TSD code for transonic unsteady aerodynamic and 
aeroelastic analysis of realistic aircraft configurations. The present capability has the option of 
half-span modeling for symmetric cases or full-span modeling to allow the treatment of 
antisymmetric mode shapes, fuselage yaw, or unsymmetric configurations such as an oblique wing or 
unsymmetric wing stores. In the present coding of the AF algorithm, the time derivatives are 
implemented for variable time stepping to allow for step-size cycling to accelerate convergence to 
steady state. Also, since the L 4 , Ln, and L; operators only contain derivatives in their respective 
coordinate directions, all three sweeps of the solution procedure have been fully vectorized. 


Results and Discussion 

Results are presented for several configurations to demonstrate and evaluate the modifications to 
the AF algorithm of the CAP-TSD code. Calculations are also presented from a flutter analysis of a 
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45° sweptback wing to assess the aeroelastic capability of the code. 


Flat Plate Airfoil Results 

Unsteady results were obtained for a flat plate airfoil at M = 0.85 to test the nonreflecting far 
field boundary conditions. 6 The flat plate airfoil was selected to allow direct comparison of results 
with the exact kernel function method of Bland. 26 The boundary conditions were tested by computing 
the lift coefficient due to the airfoil pitching about the quarter chord. Such unsteady forces are 
typically determined by calculating several cycles of forced harmonic oscillation with the last cycle 
providing the estimate of the forces. Alternatively, the forces may be obtained indirectly from the 
response due to a smoothly varying exponentially shaped pulse. 27 In this procedure, the airfoil is 
given a small prescribed pulse in a given mode of motion (in this case pitching) and the aerodynamic 
transients calculated. The harmonic response is obtained by a transfer-function analysis using fast 
Fourier transforms. Use of the pulse transfer-function technique gives considerable detail in the 
frequency domain with a significant reduction in cost over the alternative method of calculating 
multiple oscillatory responses. For the flat plate airfoil, pulse transient calculations were 
performed using 1024 time steps with At = 0.2454. The amplitude of the pulse was 0.5°. The grid 
extended 25 chordlengths above and below the airfoil, and 20 chordlengths upstream and downstream 
of the airfoil. Parallel results were obtained using reflecting (steady-state) and nonreflecting far 
field boundary conditions as shown in Fig. 1. The results are plotted as real and imaginary 
components of the unsteady lift-curve slope c, as a function of reduced frequency k. Computations 

using the reflecting boundary conditions, shown in Fig. 1(a) produce oscillations in both the real 
and imaginary parts for 0 < k < 0.2. The oscillations are produced by reflected disturbances which 
propagate back into the near field and contaminate the solution. When the calculation was repeated 
using the nonreflecting boundary conditions, shown in Fig. 1(b), the oscillations no longer occur 
since the boundary conditions absorb most of the disturbances that are incident on the grid 
boundaries. Furthermore, these results are in excellent agreement with calculations from the 
kernel function method of Ref. 26. 


F-5 Wing Results 

Calculations were next performed for the F-5 wing 28 to assess the algorithm modifications 6 to 
CAP-TSD. The F-5 wing has an aspect ratio of 3.16, a leading edge sweep angle of 31.9°, and a taper 
ratio of 0.28. The airfoil section of the wing is a modified NACA 65A004.8 airfoil which has a 
drooped nose and is symmetric aft of 40% chord. The F-5 calculations were performed using a 
constant step size for a total of 500 steps. The freestream Mach number was selected as 0.9 and the 
wing was at 0° angle of attack. The results were obtained to study the steady-state convergence 
characteristics of the modified AF algorithm. The results are presented in the form of convergence 
histories and the number of supersonic (NSUP) points versus the iteration number. 

In the original AF algorithm of Ref. 5, the Murman type-dependent switch was used. Results 
obtained using the unmodified code are presented in Fig. 2. The steady-state convergence is shown in 
Fig. 2(a); the number of supersonic points (NSUP) normalized by the final value are shown in Fig. 
2(b). For aeroelastic analysis where airloads are required rather than pressures, the solution is 
considered to be converged to engineering accuracy when a three to four order-of-magnitude 
reduction in the solution error is obtained. The "error" in the convergence history, as defined 
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herein, is the ratio of the maximum |A<|>| after n iterations to the maximum |A<|>| in the initial 
solution (first iteration). Two sets of results are plotted corresponding to two values of step size, 
At - 0.1 and 0.5. For At = 0.1 , the rate of convergence is slow and the number of supersonic points 
oscillates about the final value. Increasing the step size to At = 0.5 improves the rate of convergence 
and the oscillations in NSUP are significantly damped. The results for At = 0.5 also indicate that the 
number of supersonic points is initially more than four and one-half times the final value and that 
"spikes" begin to appear in the convergence history after 150 steps. These spikes, which represent 
a numerical instability, are due to a large transient caused by the impulsive start from a uniform 
stream using a large step size. If the calculations were started with a smaller step size, and then the 
step size increased to the larger value, the numerical instability can be avoided. Also, as shown in 
Refs. 2 and 5, the step size may be cycled through very large values such as At = 5.0 to achieve 
faster convergence to steady-state. 

The F-5 calculations with At = 0.5 were then repeated with the E-0 switch replacing the 
Murman switch. These results are labeled "unsteady algorithm" in Fig. 3. The curves are almost 
identical to the At = 0.5 curves of Fig. 2 except that the spikes in the convergence history are absent. 
The E-0 switch is more robust than the Murman switch and thus the calculation remains stable. 
Furthermore, the rate of convergence to steady-state could be increased by deleting the time 
derivatives in the residual, over-relaxing the residual, or deleting all of the time derivatives in the 
TSD equation.6 By deleting all of the time derivatives, the algorithm solves the steady TSD equation 
and is therefore referred to as the "steady algorithm." Results obtained using this algorithm are 
compared with the unsteady algorithm results in Fig. 3. The convergence history computed using the 
steady algorithm is monotonically decreasing and very smooth in comparison with the unsteady 
algorithm convergence history. The steady algorithm solution converges faster and does not produce 
the large initial overprediction of NSUP that is characteristic of the unsteady algorithm. The 
number of supersonic points converges rapidly to within 2% of its final value in only approximately 
25 steps. Over-relaxing the residual of the steady algorithm also further accelerates the 
convergence to steady-state (not shown). 


ONERA M6 Wing Results 

To test the accuracy of the modified CAP-TSD algorithm, 6 calculations were performed for the 
ONERA M6 wing.29 The M6 wing has an aspect ratio of 3.8, a leading edge sweep angle of 30°, and a 
taper ratio of 0.562. The airfoil section of the wing is the ONERA "D" airfoil which is a 10% 
maximum thickness-to-chord ratio conventional section. The freestream Mach number was selected 
as M = 0.84 and the wing was at 3.06° angle of attack. These conditions were chosen for comparison 
with the tabulated experimental pressure data of Ref. 29. This rather well-known case is a very 
challenging one, especially for a TSD code, because of the complex double shock wave which occurs on 
the upper surface of the wing. 

Steady-state calculations were performed for the M6 wing by using the AF algorithm with the 
E-0 switch. The results were obtained by cycling the step size through values as large as At = 2.0 
for a total of 500 steps. This relatively large step size corresponds to two root chords of travel per 
time step. A comparison of the resulting CAP-TSD pressures with the experimental pressure data is 
given in Fig. 4 for fj = 0.44. The data indicate that there is a relatively weak highly swept 
supersonic-to-supersonic shock wave which forms forward near the leading edge. The primary 
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supersonic-to-subsonic shock which occurs in the midchord region of the wing, coalesces with the 
first shock in the spanwise direction. Outboard toward the tip, the two shocks merge to form a single 
supersonic-to-subsonic shock wave. The CAP-TSD results obtained using first-order-accurate 
differencing in supersonic regions are in fairly good agreement with the data in predicting the 
overall pressure levels, although differences occur in the regions of the shocks. In general, the 
leading-edge suction peak is well predicted but the supersonic-to-supersonic shock is smeared. 
When the calculation was repeated using the second-order-accurate spatial differencing, a 
significant improvement was obtained in the accuracy of the results. The comparisons in Fig. 4 show 
that the supersonic-to-supersonic shock is much more sharply captured by the second-order 
method and consequently the calculated pressures are now in very good agreement with the 
experimental data. Calculations were also performed for the M6 wing using the original algorithm 
with the Murman switch. These calculations were unsuccessful because of a numerical instability 
which was produced by the highly expanded flow about the leading edge of the wing. 

An unsteady calculation was also performed for the M6 wing at M = 0.84, to investigate the 
robustness of the modified algorithm for time-dependent applications. In this demonstration 
calculation, the wing was forced to oscillate in pitch about a line perpendicular to the root at the root 
midchord. The amplitude of the motion was 2° peak-to-peak about the mean angle of attack of ao = 
3.06°. The reduced frequency was selected as k = 0.1 and only 300 steps per cycle of motion were 
used. This corresponds to a step size of At = 0.1047. Three cycles of motion were computed to 
obtain a periodic solution. Unsteady pressure distributions, obtained using first-order and second- 
order accurate supersonic differencing, are shown in Fig. 5 at the maximum pitch angle (a « 4.06°) 
for n = 0.44. Similar to the steady-state results, these pressure comparisons illustrate that the 
supersonic-to-supersonic shock is more sharply captured by the second-order method. Further 
instantaneous pressure distributions at two points during the third cycle of motion are shown in Fig. 
6 for five span stations along the wing. Pressures at the wing maximum angle of attack (a = 4.06°) 
and pressures at the wing minimum angle of attack (a = 2.06°) are both presented in the figure. As 
the wing pitches up, the shocks move aft and the supersonic-to-subsonic shock grows in strength. 
As the wing pitches down, the shocks move forward and the supersonic-to-supersonic shock is more 
sharply defined. For this case, both of the shocks oscillate over approximately 10% of the chord 
during a cycle of motion. Also, the supersonic-to-supersonic shock at fi = 0.80 periodically appears 
and disappears during a cycle of motion. The results illustrate the large shock motions that the 
modified AF algorithm is capable of computing. The improved algorithm captures the shocks sharply 
and is sufficiently robust to compute this complex unsteady flow using only 300 steps per cycle of 
motion. 

To test the entropy and vorticity modifications to TSD theory ,7 further calculations were 
performed for the M6 wing. Results were obtained at the freestream Mach number of M = 0.92 with 
the wing at 0° mean angle of attack. These conditions correspond to an AGARD test case for 
assessment of inviscid flow field methods 30 and were selected for comparison with the tabulated 
Euler results of Rizzi contained therein. Calculations were performed using: (a) unmodified TSD 
theory and (b) TSD with entropy and vorticity effects. Steady pressure distributions along three 
span stations (rj = 0.08, 0.47, and 0.82) of the wing are presented in Fig. 7 from both solutions. 
For this steady-state case, the flow is symmetric about the wing with shocks on the upper and lower 
surfaces. As shown in Fig. 7(a), the results from the unmodified TSD theory compare well with the 
Euler results in predicting the leading-edge suction peak and the overall pressure levels. However, 
the shock is located too far aft and is too strong outboard near the tip (at rj = 0.82, for example), in 
comparison with the Euler calculation. When the entropy and vorticity modifications are included in 
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the calculation, the shock is displaced forward from the previous solution, as shown in Fig. 7(b). 
Here the shock location and shock strength are in very good agreement with the Euler results at all 
three span stations on the wing. Consequently, the steady pressure distributions from the modified 
TSD theory now compare very well with the Euler pressures. 


General D ynamics F- 16 C Aircraft Model Results 

Results were also obtained for the General Dynamics F-16C aircraft model 31 to investigate 
application of the modified algorithm to a realistic aircraft configuration. 3 Shown in Fig. 8 are the 
F-16C components that are modeled using CAP-TSD. The F-16C is modeled using four lifting 
surfaces and two bodies. The lifting surfaces include: (1) the wing with leading and trailing edge 
control surfaces, (2) the launcher, (3) a highly swept strake, aft strake, and shelf surface, and 
(4) the horizontal tail. The bodies include: (1) the tip missile, and (2) the fuselage. In these 
calculations, the freestream Mach number was M = 0.9 and the F-16C aircraft was at 2.38° angle of 
attack. Also, the leading-edge control surface of the wing was deflected upwards 2° for comparison 
with the experimental steady pressure data of Ref. 32. Furthermore, the calculations were 
performed on a grid which conforms to the leading and trailing edges of the lifting surfaces and 
contains 324,000 points. Since the grid is Cartesian, it was relatively easy to generate, even for 
such a complex configuration as the F-16C aircraft. Also, the calculations required only about 0.88 
CPU seconds per time step and thirteen million words of memory on the CDC VPS-32 computer at 
NASA Langley Research Center. 

Steady-state calculations were performed for the F-16C aircraft using the AF algorithm with the 
E-0 and Murman switches. The E-0 results were obtained using both the first-order and second- 
order accurate supersonic differencing. Steady pressure comparisons are given in Fig. 9 for three 
span stations of the wing and one span station of the tail. Both sets of E-0 results are presented for 
comparison with the experimental data. The results obtained using the Murman switch were 
originally published in Ref. 2. These results are identical to plotting accuracy with the first-order 
E-0 results, and therefore are not shown. The steady pressure comparisons indicate that there is a 
moderately strong shock wave on the upper surface of the wing and the CAP-TSD pressures agree 
well with the experimental pressures. For the tail, the flow is predominantly subcritical and the 
calculated results again agree well with the data. Comparison of pressures computed using first- 
order and second-order accurate supersonic differencing shows very small differences. The largest 
difference, for example, occurs on the wing at iiw = 0.79 where the second-order calculation 
predicts a slightly stronger shock. 

Unsteady results were also obtained for the F-16C aircraft to investigate the robustness of the 
modified algorithm for realistic-aircraft time-dependent applications. For simplicity, the 
calculation was performed for a rigid pitching motion where the entire aircraft was forced to 
oscillate about the model moment reference axis at a reduced frequency of k = 0.1. The oscillation 
amplitude was chosen as ai = 1.5° which is three times the value used to obtain similar results 
presented in Ref. 2. Three cycles of motion were computed using 300 steps per cycle of motion 
corresponding to At = 0.1047. Calculations were performed using both the Murman and E-0 
switches. The solution using the original algorithm with the Murman switch, however, was 
numerically unstable for this case as shown in Fig. 10. Instantaneous pressure distributions at time 
steps 94 and 95 are plotted in the figure, computed using the Murman (Fig. 10(a)) and E-0 (Fig. 
10(b)) switches. The numerical instability begins in the region of the launcher/tip-missile where 
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the grid spacing is smallest. Figure 10(a) shows the instability in the form of an oscillation in the 
wing upper surface pressure distribution at ijw = 0.94 from approximately 30% to 60% chord. 
The program subsequently failed during step 96 which is 21 steps after the maximum pitch angle in 
the first cycle of motion. The calculation involving the modified algorithm (E-0 switch with the 
first-order accurate supersonic differencing) is stable, however, as shown in Fig. 10(b). Here the 
pressure distributions for steps 94 and 95 are very similar and the calculation proceeds with no 
difficulty. In fact, the modified AF algorithm with the E-0 switch is numerically stable for this case 
with either the first-order or second-order supersonic differencing. 

Unsteady pressure distributions along the wing and tail during the third cycle of motion are 
shown in Fig. 11, computed using the E-0 switch with the second-order accurate supersonic 
differencing. Two sets of calculated pressures are presented corresponding to the aircraft at the 
maximum (a = 3.88°) and minimum (a = 0.88°) pitch angles. Comparison of the results indicates 
that large changes in pressure occur along the upper and lower surfaces of the wing as the aircraft 
oscillates in pitch. For example, the shock on the wing upper surface oscillates over more than 10% 
of the chord during a cycle of motion. Also, the shock is approximately twice as strong at the 
maximum pitch angle as it is at the minimum pitch angle. For the tail, the changes in the pressure 
distributions due to aircraft pitching are relatively very small in comparison with the changes in 
wing pressures, as further shown in Fig. 11. The tail is located considerably aft of the pitch axis and 
thus its motion is plunge dominated which results in much smaller airloads for the low value of k 
considered. 


Wing Flutter Results 

To assess the CAP-TSD code for flutter applications, a simple well-defined wing case was selected 
as a first step toward performing aeroelastic analyses for complete aircraft configurations .4 The 
wing being analyzed is a semispan wind-tunnel-wall-mounted model that has a quarter-chord sweep 
angle of 45°, a panel aspect ratio of 1 .65, and a taper ratio of 0.66. 33 The wing is a proposed 
AGARD standard aeroelastic configuration which was tested in the Transonic Dynamics Tunnel (TDT) 
at NASA Langley Research Center. 34 A planview of the wing is shown in Fig. 12. The wing has a 
NACA 65A004 airfoil section and was constructed of laminated mahogany. In order to obtain flutter 
for a wide range of Mach number and density conditions in the TDT, holes were drilled through the 
wing to reduce its stiffness. To maintain the aerodynamic shape of the wing, the holes were filled 
with a rigid foam plastic. A photograph of the model mounted in the TDT is shown in Fig. 13. The 
wing is being modeled structurally using the first four natural vibration modes which are 
illustrated in Figs. 14 and 15. Figure 14 shows oblique projections of the natural modes while Fig. 
15 shows the corresponding deflection contours. These modes which are numbered 1 through 4 
represent first bending, first torsion, second bending, and second torsion, respectively, as 
determined by a finite element analysis. The modes have natural frequencies which range from 9.6 
Hz for the first bending mode to 91.54 Hz for the second torsion mode. 

Flutter calculations were performed for the 45° sweptback wing using CAP-TSD to assess the 
code for aeroelastic applications. Two sets of results are presented corresponding to: (1) using the 
linear potential equation (F = G = H = 0) and modeling the wing aerodynamically as a flat plate (zero 
thickness) and (2) using the complete (nonlinear) TSD equation and including wing thickness. The 
first set of results allows for direct comparison with parallel linear theory calculations performed 
using the FAST subsonic kernel function program. 35 The second set of results more accurately 
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models the wing geometry as well as the flow physics. All of the results are compared with the 
experimental flutter data of Ref. 34 which spans the range, 0.338 < M < 1.141. 

Comparisons of flutter characteristics from the linear calculations with the experimental data 
are given in Fig. 16. Plots of flutter speed index (defined as U/(bo coa Vji)) and nondimensional 

flutter frequency (defined as ©/wa) as functions of freestream Mach number, are shown in Figs. 

16(a) and 16(b), respectively. The experimental flutter data defines a typical transonic flutter 
"dip" with the bottom near M = 1.0 for this case. (Note that these results are shown with an 
expanded scale.) The bottom of the dip in flutter speed index (Fig. 16(a)) was defined by the 
approach to the M = 1.072 flutter point during the wind tunnel operation. Results from the CAP- 
TSD (linear) code are presented at twelve values of M covering the entire Mach number range over 
which the flutter data was measured. Results from the FAST program are presented for the limited 
range 0.338 0.96 since the method is restricted to subsonic freestreams. Overall, the linear 

CAP-TSD results compare well with the experimental data for subsonic as well as supersonic Mach 
numbers. Note that the subsonic FAST results are also in good agreement with the data. Such a result 
is not unexpected for this very thin wing of moderate sweep and taper at zero angle of attack. It does 
indicate that the wing properties are well-defined for benchmark purposes. 

In the subsonic Mach number range, the CAP-TSD and FAST calculations predict a slightly 
unconservative flutter speed, except at M = 0.338, by as much as 2% (Fig. 16(a)), and a higher 
flutter frequency (Fig. 16(b)) in comparison with the experimental data. In general though, the 
linear CAP-TSD results agree well with the FAST results in both flutter speed and frequency. The 
good agreement in this three-way correlation between experiment, linear theory, and CFD flutter 
results gives confidence in the CAP-TSD code for flutter prediction. 

Comparisons of flutter characteristics from the linear and nonlinear CAP-TSD calculations with 
the experimental data are given in Fig. 17. Figure 17(a) shows flutter speed index versus Mach 
number and Fig. 17(b) shows nondimensional flutter frequency versus Mach number. Three flutter 
points are plotted from the nonlinear CAP-TSD calculations corresponding to M - 0.678, 0.901 , and 
0.96. Comparisons between the two sets of CAP-TSD results show differences due to wing thickness 
and nonlinear effects. With increasing Mach number these differences become larger. For example, 
at M = 0.678, 0.901, and 0.96, the flutter speed index decreased by 1%, 5%, and 19%, 
respectively, as shown in Fig. 17(a). Similar decreases also occur in the flutter frequency (Fig. 
17(b)). The decrease in flutter speed at M = 0.901 is largely due to including wing thickness since 
there are no supersonic points in the flow at this condition. The decrease in flutter speed at M = 
0.96 is attributed to both wing thickness and nonlinear effects since an embedded supersonic region 
of moderate size was detected in the wing tip region. The nonlinear CAP-TSD results at both M = 
0.901 and 0.96 are slightly conservative in comparison with the experimental flutter speed index 
values. Nonetheless, the nonlinear CAP-TSD flutter results compare favorably with the 
experimental data, which is the first step toward validating the code for general transonic 
aeroelastic applications. 


Concluding Remarks 

A transonic unsteady aerodynamic and aeroelasticity code called CAP-TSD has been developed for 
application to realistic aircraft configurations. The name CAP-TSD is an acronym for Computational 
Aeroelasticity Erogram - Transonic Small disturbance. The code permits the calculation of steady 
and unsteady flows about complete aircraft configurations for aeroelastic analysis in the flutter 
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critical transonic speed range. The CAP-TSD code uses a time-accurate approximate factorization 
(AF) algorithm for solution of the unsteady transonic small-disturbance equation. An overview of 
the code development effort was given and recent algorithm modifications were described. The 
algorithm modifications include: an Engquist-Osher (E-O) type-dependent switch to treat regions of 
supersonic flow, extension of the E-0 switch for second-order spatial accuracy, nonisentropic 
effects to treat strong-shock cases, nonreflecting far field boundary conditions for unsteady 
applications, and several modifications to accelerate convergence to steady-state. Calculations were 
presented for several configurations including the General Dynamics one-ninth scale F-16C aircraft 
model to evaluate the algorithm modifications. The modifications have significantly improved the 
stability of the AF algorithm and the reliability of the CAP-TSD code in general. 

Results were also presented from a flutter analysis of a 45° sweptback wing. The flutter 
boundaries from CAP-TSD (linear) were in agreement with parallel subsonic linear theory results 
and compared well with the experimental flutter data for subsonic and supersonic freestream Mach 
numbers. The preliminary nonlinear CAP-TSD flutter results also compared favorably with the 
experimental data which is the first step toward validating the code for general transonic aeroelastic 
applications. 
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Reduced frequency k 


(a) reflecting boundary conditions. 



Reduced frequency k 


( b ) nonreflecting boundary conditions. 

Fig. 1 Comparisons of unsteady lift-curve slope for a flat plate airfoil at M = 0.85. 
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( a ) steady-state convergence. 
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( b ) number of supersonic points. 


Effects of step size on the solution computed using the Murman switch for the F-5 
wing at M = 0.9 and ao = 0°. 




( a ) steady-state convergence. 



( b ) number of supersonic points. 

Fig. 3 Effects of deleting all TSD time derivatives on the solution computed using the 
Engquist-Osher switch for the F-5 wing at M = 0.9 and ao = 0°. 
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Fig. 4 Effects of first-order and second-order accurate supersonic differencing on the steady 
pressure distributions of the ONERA M6 wing at M = 0.84 and ao = 3.06°. 



X/C 

Fig. 5 Effects of first-order and second-order accurate supersonic differencing on the 
unsteady pressure distributions of the ONERA M6 wing during the third cycle of rigid 
pitching at M = 0.84, ao = 3.06°, ai = 1 .0°, and k = 0.1 . 
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Fig. 6 Instantaneous pressure distributions on the ONERA M6 wing during the third cycle of 
rigid pitching at M = 0.84, ao = 3.06°, ai = 1.0°, and k = 0.1. 



Fig. 7 Comparison of steady pressure distributions for the ONERA M6 wing at M = 0.92 and 
ao = 0°. 
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+ vorticity 
o Euler 



(b) TSD + entropy + vorticity. 
Fig. 7 Concluded. 



Fig. 8 CAP-TSD modeling of the General Dynamics one-ninth scale F-16C aircraft model. 
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Fig. 9 Comparisons between CAP-TSD steady pressure distributions computed using first- 
order and second-order accurate supersonic differencing with the experimental 
pressure data for the wing and tail of the F-16C aircraft model at M = 0.9 and 
ao = 2.38°. 



(a) numerical instability with Murman switch. 

Fig. 10 Effect of type-dependent switch on numerical stability for rigid pitching of the F-16C 
aircraft model at M = 0.9, ao = 2.38°, ai = 1.5°, and k = 0.1. 
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(b) improved numerical stability with Engquist-Osher switch. 


Fig. 10 Concluded. 



Fig. 11 Instantaneous pressure distributions on the wing and tail of the F-16C aircraft model 
during the third cycle of rigid pitching at M = 0.9, ao = 2.38°, ai = 1.5°, and 
k = 0.1. 
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Fig. 12 Planview of 45° sweptback wing. 



Fig. 13 45° sweptback wing in the NASA Langley Transonic Dynamics Tunnel. 
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( a ) flutter speed index versus Mach number. 
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( b ) nondimensional flutter frequency versus Mach number. 


Fig. 16 Comparisons of linear flutter calculations with experimental data for the 45° 
sweptback wing. 
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(a) flutter speed index versus Mach number. 
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( b ) nondimensional flutter frequency versus Mach number. 


Fig. 17 Comparisons of linear and nonlinear CAP-TSD flutter predictions with experimental 
data for the 45° sweptback wing. 
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